1. Field of the Invention
The present invention relates to a device simulation method for use in numerical analyses of a semiconductor device and, more particularly, to a device simulation method effective in a case where some area of a semiconductor device as an object of processing includes a part at which drastic change of an electric field occurs.
2. Description of the Related Art
In conventional numerical analyses of semiconductor devices, drift-diffusion models are widely used in which carriers (electron and hole) are approximated as fluid. Although energy transport models with a high level of approximate are used in some cases, description will be here made of the drift-diffusion models for the convenience of explanation.
In device simulation using a drift-diffusion model in the steady-state, such equation of conservation of charge, electronic current continuity equation and hole current continuity equation as shown below are set as basic equations. These basic equations are detailed, for example, in the literature "Process Device Simulation Technique" (written and edited by Ryo Dan, pp. 99-105, Sangyo Tosho, 1990). EQU div DV (V denotes vector)=.rho.(equation of conservation of charge)(1) EQU DV=.epsilon.EV (2) EQU EV=-grad .psi. (3) EQU .rho.=q(p-n+N.sub.D -N.sub.A) (4)
where D denotes dielectric flux density, .rho. charge density, EV electric field, .epsilon. dielectric constant, q elementary charge, p hole density, n electron density, N.sub.D donor density and N.sub.A acceptor density. EQU div JV.sub.n =q(R-G)(electron current continuity equation) (5) EQU div JV.sub.p =-q(R-G)(hole current continuity equation) (6)
where JV.sub.n denotes electron current, JV.sub.p hole current, R carrier re-combination term and G carrier generation term. EQU JV.sub.n =q.multidot.n.multidot..mu..sub.n .multidot.EV+q.multidot.D.sub.n .multidot.grad n (7) EQU JV.sub.p =q.multidot.n.multidot..mu..sub.p .multidot.EV-q.multidot.D.sub.p .multidot.grad p (8)
where .mu..sub.n denotes electron mobility, .mu..sub.p hole mobility, D.sub.n electron diffusion coefficient and D.sub.p hole diffusion coefficient. EQU D.sub.n =.mu..sub.n .multidot.(k.sub.B .multidot.T/q) (9) EQU D.sub.P =.mu..sub.p .multidot.(k.sub.B .multidot.T/q) (10)
where k.sub.B denotes Boltzmann's constant and T denotes temperature.
In the above equations, variables to be solved are potential .psi., electron density n and hole density p. In general, with a plurality of designated biases as boundary conditions, the biases are sequentially updated to calculate the equation of conservation of charge (1), the electron current continuity equation (5) and the hole current continuity equation (6).
In actual numerical calculation, these basic equations are discretized and calculated on predetermined analysis meshes.
One of essential performance barometers of a semiconductor device is a breakdown voltage. When a high voltage is applied to a semiconductor device, the strong electric field generated in the semiconductor causes conduction electrons to gain high energy to collide with electrons in a valence band, whereby the colliding electrons are excited by a conductor to further collide with electrons in other valence band. Repetition of such collision and excitation results in generation of a large quantity of pairs of electrons and holes by geometrical progression to rapidly increase current. An applied voltage involving such phenomenon is called a breakdown voltage. In addition, the rapid generation of electrons and holes is called avalanche amplification or impact ionization.
As a breakdown voltage estimating manner by numerical calculation, there is a method of obtaining a multiplication factor of current by conducting ionization integration along a line of electric force. This technique is recited, for example, in the literature "Physics of Semiconductor Devices (Second Edition)" (written by S. M. Sze, pp. 99, John Willy & Sons Inc., 1981). The ionization integration in this literature is represented by the right side of the following expression and the multiplication factor is a value represented by M.sub.p of the left side of the following expression. ##EQU1## where x and x' represent integration variables of a one-dimensional coordinate along a line of electric force, W represents a range of integration denoting a coordinate at the end of the line of electric force, and .alpha..sub.P and .alpha..sub.n represent impact ionization coefficients of hole and electron.
An applied bias making the multiplication factor M.sub.P infinite is referred to as a breakdown voltage. The multiplication factor goes infinite in the equation (11) under the condition that the ionization integration on the right side of the equation takes the integration value of 1.
In addition, various models are proposed for the impact ionization coefficients .alpha..sub.P and .alpha..sub.n. For example, in the well-known Sze and Clowell model, the coefficients are expressed by the following equations (12) to (15) as recited in the literature "Temperature Dependence of Avalanche Multiplication in Semiconductors" (C. R. Crowell, S. M. Sze, Applied Physics Letters, vol. 9, No. 6, pp. 242-244, 1966). Parameters of the equations (12) to (15) are shown in Table 1. ##EQU2## where E represents electric field strength and T represents temperature.
TABLE 1 ______________________________________ (electron) (hole) parameter value unit Parameter value unit ______________________________________ a.sub.e 0.426 1/V! a.sub.h 0.243 1/V! b.sub.e 4.81E5 V/cm! b.sub.h 6.53E5 V/cm! c.sub.e 3.05E-4 1/K! c.sub.h 5.35E-4 1/K! d.sub.e 6.86E-4 1/K! d.sub.h 5.67E-4 1/K! ______________________________________
For conducting the above-described numerical calculation along a line of electric force, the line of electric force should be obtained in advance. More specifically, an electric field distribution is calculated based on a potential distribution obtained by solving the basic equations and a line of electric force is traced with respect to the obtained electric field distribution to obtain a desired line of electric force. The electric field then is a vector quantity having a value on a nodal point of a mesh used for solving the basic equations. Conventional device simulation methods for obtaining lines of electric force are disclosed, for example, in the literature "Personal Computer Electromagnetics" (Takehiko Adachi, Hitoshi Sekimoto, Yasuaki Watanabe, Asakura Books, 1986, pp. 30-31) and the literature "Simulation of Electromagnetic Phenomenon Learned by Personal Computer" (Noriyuki Kitahara, Shinnichi Hirata, Morikita Publishing, 1991, pp. 12-13). The conventional device simulation methods recited in these literatures simply interpolate and trace values of electric fields based on a value of an electric field on each nodal point of the mesh.
FIG. 16 is a flow chart showing a conventional device simulation method. With reference to FIG. 16, the conventional device simulation method includes Step 1603 for conducting interpolation based on electric fields at both end points of a mesh edge where the point "m.sub.i " (i=1, 2, 3, . . . ) is located to obtain a search direction vector "dmV.sub.i " (hereinafter denoted as "dm.sub.i ") of the line of electric force, Step 1604 for extending a line segment from the point "m.sub.i " in the direction of the search direction vector "dm.sub.i " of the line of electric force to obtain a point "m.sub.i+1 " of intersection between the line segment and a mesh edge, and Step 1605 for determining whether the intersection point reaches the end point of the line of electric force.
With reference to the flow chart illustrated in FIG. 16 and FIG. 17 which schematically shows a processing procedure according to the conventional device simulation method, operation of the conventional device simulation method will be described. Triangular meshes represented by dotted lines in FIG. 17 are a part of the meshes to be used for solving the basic equations. Electric field at each mesh nodal point "n.sub.i " (i=1-9) is represented by vector "E(EV)n.sub.i " (i=1-9). Under the initial conditions, a point "m.sub.1 " on the mesh edge "n.sub.1 n.sub.2 " is regarded as a starting point of a line of electric force.
First, calculate a search direction vector of the line of electric force (Step 1603). In a case where the line of electric force is traced in the forward direction, that is, in the potential gradient descending direction, a search direction vector "dm.sub.1 " of the line of electric force at the starting point "m.sub.1 " will be obtained by the following expressions through linear interpolation of the electric field "Em.sub.1 " from the vectors "En.sub.1 " and "En.sub.2 " to "m.sub.1 ". ##EQU3## where l.sub.n1n2 represents a length of a mesh edge between n.sub.1 and n.sub.2, l.sub.n1m1 represents a length of a mesh edge between n.sub.1 and m.sub.1 and l.sub.n2m1 represents a length of a mesh edge between n.sub.2 and m.sub.1.
Next, find a point "m.sub.2 " at which a line segment extended from the starting point "m.sub.1 " in the direction of the search direction vector "dm.sub.1 " of the line of electric force intersects with a mesh edge (Step 1604). In FIG. 17, this intersection point "m.sub.2 " is located on the mesh edge n.sub.2 n.sub.4. Then, determine whether the point "m.sub.2 " is the end point of the line of electric force (Step 1605). Next, find a search direction vector "dm.sub.2 " of the line of electric force at the point "m.sub.2 " and find a point "m.sub.3 " at which a line segment extended in the direction of the vector intersects with a mesh edge. Then, determine whether the point "m.sub.3 " is the end point of the line of electric force (Step 1605).
Points "m.sub.4 ", "m.sub.5 ", "m.sub.6 " and "m.sub.7 " can be found in the same manner. As a result, such a line of electric force "m.sub.1, m.sub.2, m.sub.3, m.sub.4, m.sub.5, m.sub.6, m.sub.7 " as illustrated in the figure is obtained. Every time each point is obtained, determination is made whether the point is the end point of the line electric force or not and if, for example, the point "m.sub.7 " is the end point, processing is completed. At Step 1605, the point as an object of the determination is regarded as reaching the end point of the line of electric force when a traced line of electric force reaches an outline part of the meshes or when an electric field becomes weak as the potential distribution goes flat.
When the line of electric force is traced in the reverse direction from the starting point, calculation is made by the following expression, with a sign of a search direction vector of the line of electric force changed. ##EQU4##
Description will be next made of processing to be conducted according to the conventional device simulation method when a device as an object of simulation includes an electric field drastically varying region, with reference to FIG. 18(A) and (B) which schematically show the processing procedures. In the figures, a linear part represented by a chain dotted line is a valley-formed potential part where the potential changes in the form of a valley. Such drastic change of electric field occurs, for example, in the vicinity of a drain right under a gate oxide film when a voltage is applied across a gate and a drain in a MOS transistor. Also at a region where such electric field changes as shown in FIG. 18(A) occurs, like the case described with reference to FIG. 17, a search direction vector "dm.sub.1 " of the line of electric force is first calculated regarding a point "m.sub.1 " as a starting point of a line of electric force to find a point "m.sub.2 " at which a line segment extended in the direction of the vector "dm.sub.1 " intersects with a mesh edge. Next, find a search direction vector "dm.sub.2 " of the line of electric force at the point "m.sub.2 " to find a point "m.sub.3 " at which the line segment extended in the direction of the vector intersects with a mesh edge. By tracing "m.sub.4 ", "m.sub.5 ", "m.sub.6 ", "m.sub.7 ", "m.sub.8 " and "m.sub.9 " in the same manner, a line of electric force m.sub.1, m.sub.2, m.sub.3, m.sub.4, m.sub.5, m.sub.6, m.sub.7, m.sub.8 and m.sub.9 is obtained.
In this example, the line of electric force is obtained as a line fluctuating with the valley-formed potential part in between. Actual line of electric force, however, extends along the valley-formed potential part after reaching the point "m.sub.2 " of intersection with the valley-formed potential part, as represented by the solid line in FIG. 18(B). The line of electric force to be used in numerical calculation will be m.sub.1, m.sub.2, m.sub.3 represented by the solid line. When numerical calculation is made along such a line of electric force fluctuating as shown in FIG. 18(A) and obtained by the conventional device simulation method, calculation error would be therefore increased.
The reason here for using in numerical calculation not the line of electric force represented by the two dot chain line which is an actual line of electric force but the line of electric force represented by the solid line which is an approximate line of electric force is that the latter makes calculation speed higher and data handling in a program easier than the former does. Moreover, since a unit of discretization is specified by a mesh generated at a semiconductor device, such approximation as represented by the solid line in FIG. 18(B) will involve no serious calculation error in ionization integration calculation.
As described in the foregoing, the conventional device simulation method has a shortcoming that when a device to be simulated includes an electric field drastically changing portion such as a valley-formed potential part, a search direction vector of a line of electric force fluctuates, so that an incorrect line of electric force is recognized to increase calculation error.